An analytically solvable model of probabilistic network dynamics 
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We present a simple model of network dynamics that can be solved analytically for uniform 
networks. We obtain the dynamics of response of the system to perturbations. The analytical 
solution is an excellent approximation for random networks. A comparison with the scale-free 
network, though qualitatively similar, shows the effect of distinct topology. 
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Recent advances in the understandingof complex so- 
cial 0, biological 0, and technological |3] systems have 
revealed widespread if not universal properties of the 
topology of networks of association, interaction and com- 
munication. These properties, include small-world global 
connectivity , scale-free local connectivity distribution 
0, and characteristic local motif structures!^] . Central 
to our understanding of complex systems Q is charac- 
terizing their response to environmental stimuli. While 
much of the focus has been on robustness to random per- 
turbation or directed attack, Q the effectiveness of re- 
sponse requires satisfying a wider range of requirements 
including, for example, sensitivity to particular stimuli 
Indeed, one of the main functions of biological and 
social systems is the detection of specific stimuli that 
require collective (large scale) response in seeking desir- 
able resources (foraging) or responding to dangers ( "fight 
or flight"). Thus understanding the nature of system 
function and behavior from topological structure requires 
mapping the interaction structure given by a topologi- 
cal network onto the dynamics of system response [lOj. 
Insofar as the network of interactions obtained in re- 
cent research characterizes the internal interactions of 
a system, these interactions must provide key informa- 
tion about the dynamics of response to external pertur- 
bations. Therefore, the construction of solvable dynamic 
models is essential for understanding the general features 
of the problem. In this Letter we propose such a model 
of probabilistic network dynamics, and we solve it ana- 
lytically for uniform networks. 

We consider a general network with N nodes. To each 
node i is assigned an internal state <ii that can take the 
values or 1. At each time step the state of a node is 
updated according to the following rule: either the state 
does not change, which happens with probability p or, 
with probability (1 — p), it copies the state of one of 
its neighbors. This process describes, for example, the 
behavior of a group of high school students choosing to 
adopt one style of dress or another, or the propagation of 
a mutation through a species ^lj. Since the states of a 
node are abstract labels, the change of one node to adopt 
the state of another can be considered a general model 
of influence propagation, with each node state a label for 
its own relevant physical property. 



The 2 N states of the network can be labeled by a string 
of zeroes and ones describing the internal state of each 
node in sequence (ctn-io-n-2 ■ ■ ■o~x<jq). Alternatively, the 
states can be labeled by integers via x — Ylj=o &j2? , 
with x varying between and 2 N — 1. 

Let Pt{x) be the probability of finding the network in 
the state x at time t and let the network evolve through 
asynchronous updates, where a single node is allowed to 
change at each time step. To find how this probability 
changes with time we define the auxiliary state Xk which 
is equal to x at all nodes except at node k, which has the 
opposite internal state. The probability of finding the 
network in the state x at time t + 1 can now be written 
as a sum of three terms: (a) the probability that the 
network was in state x at time t and that the selected 
node did not change plus (b) the probability that it was 
in the state x and the selected node copied the state of 
an identical neighbor plus (c) the probability that the 
network was in the state Xk at time t and that the node 
k was selected and its state a k — 1 — c/c changed to a k ■ 



P t+ i(x) =pP t {x) + 



N 



^2{P t (x) Prob[cj k -^<J k } 



+ P t (x k ) Prob[a k > } • 

The probability Prob[<7k — > erj is just the number of 
neighbors of node k in the state a k divided by the total 
number of neighbors (the degree) d k = J^iLo 1 ^ik, where 
Cik is the connectivity (or adjacency) matrix. This can 
be written as 



JV-l 



dk 



i=0 



The probability Prob[<7k — * Ok] is also given by this for- 
mula, since a k = 1 — a k and d~i = Gi for i ^ k. Using 
these relations, we obtain the following master equation 
for the network dynamics: 



P t+ i{x) = pP t (x) + 



N 



N-l 



N-l 



k=0 i=0 



(1) 



\l-<Ti-a k \[Pt{x)+Pt(x k )] 



Finding Pt(x) for networks with arbitrary topologies 
can be very difficult. However, the problem can be 



2 



completely solved for fully connected networks, where 
dk = N — 1. In this case the nodes are indistinguish- 
able from each other and the states of the network can 
be labeled simply by counting the number of nodes in the 
internal state 1, given by n(x) = X^ *- The probability 
of finding the network in the state labeled by n is related 
to P(x) by 



P{n{x)) = P{x)B(N,n) 



(2) 



where B(N, n) = N\/[n\(N—ri)\] is a binomial coefficient. 
We now simplify the last two terms on the right of Eq.JIJ. 
To do this we separate the sum over k into the cases 
(jfc = 1 and Ufc = 0. For the first of these terms we obtain 



+ C ik {l - ai)]P t (x) 
(N-n){N-n- l)]P t (x) . 



For the third term we observe that P t (xk) corresponds 
to the state n — 1 if a k = 1 and to n + 1 if = 0. When 
we separate the sum over k into the cases ak — 1 and 
(Tfe = we write Pt{xk] 1) and Pt(xk\ 0) respectively. We 
obtain 



N(N-l) 



[n(n-l)P t (x k ;l)+(N-n){N-n-l)P t (x k ;0)] 



Substituting these terms into the master equation and 
multiplying both sides by B(N,n) we obtain, after some 
simplification, 



Pi 



, + i(n) =pPt(n) + ww 
{[n(n-l) + (JV 



n)(N-n- l)]Pt(n) 
+ (N - n + l)(n - l)P t (n - 1) 
+(N-n-l)(n + l)P t (n + l)} . 



(3) 



For uniform networks where = do is the same for 
all nodes, if do < N — 1 states with the same n(x) can be 
distinguished by the way the internal states with a k = 1 
are distributed among those with a k — 0. However, if 
we combine states with the same n(x), the procedure de- 
scribed above can still be applied. In this case the factor 
dk in the denominator is replaced by do. However, on 
the average (with respect to the different states labeled 
by n), the counting of the number of neighbors must be 
multiplied by do/ (JV — 1), so that do cancels and we get 
(JV — 1) back in the denominator. Therefore equation 
© holds in this case as well. For random networks the 
degree of each node is nearly constant, and we can still 
use equation Q as an approximation for the dynamics. 

We now proceed to the calculation of the transition 
probabilities. The probabilities Pt(n) define a vector P t 
of N + 1 components. The master equation can be 
written in matrix form as Pt+i = U Pt where the evo- 
lution matrix U is tridiagonal. The propagation of an 
initial probability vector requires the calculation of pow- 
ers of U . Alternatively, we can diagonalize U and use its 
eigenvectors as a basis. This approach has been used 01 
to calculate the eigenvalues of the transition matrix for 



certain population models. Here we shall calculate not 
only the eigenvalues but also the eigenvectors, obtaining 
the complete solution of the dynamical problem. 

The eigenvalues of U can be calculated for small matri- 
ces and extrapolated to matrices of arbitrary size. They 
are given by 



K = 1 - 



1-p 



N(N - 1) 



-(r - 1) . 



with r = 0,1,..., JV. The only degeneracy occurs for 
Ao = Ai = 1. The other eigenvalues are all smaller than 
1 and decrease towards Xn = V- 

Since U is not symmetric, its eigenvectors do not form 
an orthogonal set. Let \a r ) and (b r \ be the right and left 
eigenvectors of U, with components a rrn and b rm . Then 



1 



N 
r=0 



■\a r ){b r \ = 1 



where = (b r >\a r ) — T r S rr i and r r = ^ m a rm b rm . 

An initial vector |f(0)) containing the information 
about the probability of the different states at time zero 
can be projected using this resolution of unity and easily 
evolved: 



1 



\ V (t))^U t \v(0))=J2^( b r\ V (°)) X rK) 
1 r 



r=0 



The transition probability between two network states 
with n = M and n = L after a time t can now be cal- 
culated by taking the components of the initial vector as 
V m (0) = SM,m and projecting the evolved state onto the 
state with components o~L,m- 

N 1 

P(L, t; M, 0) = TT°rMarLK ■ 

r=0 r 

The coefficients a rm follow a recursion relation that 
can be derived directly from the eigenvalue equation for 
U. For r = and r = 1 the eigenvectors can be found 
immediately: 



l a o) = 



and 



IM = (in 



o 







v i / 



1 1) 



\ai) = 



1 \ 








V i / 



(4) 



|&i) = (JV JV-2 JV-4 



-7V+2 -JV) 



(5) 



In order to calculate the remaining eigenvectors, we 
define the auxiliary eigenvalues fi r by 

Pr = (1 - A r ) — — — ^j— — - — — =r(r-l) (6) 
1-p 
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and the auxiliary coefficients 

A rm = m{N — m)a rm m = 1, 2, . . . , N — 1 . (7) 

The recursion relation for the A rm can be written explic- 
itly as: 



-1 2A T ra -\- A r rn — l 



IXr ( A. 



A < ■ 



N \ m N — m 
A generating function is now defined as 



(8) 



N-l 



m=l 

(note that A r0 — A t n — 0). Multiplying eq.JEJ by x" 
and summing over to we get, on the left side, 



fr 

— (1 - xf - A r i + A r N-lX 



N 



Letting <p r (x) = a rm x m then 

4 1 



■(*)=E 



m(iV — to) 



iV 



Differentiating with respect to x, using Ea. ljf f|) . dividing 
by x N_1 and differentiating again, we find 



N-l 



1 F r 



x {1-xf 







Comparing with Eq. (|f 3|) we see that <p r — —F r /ji r . 
Therefore, except for a normalization, the generating 
function for the coefficients a rm , <p r {x), is equal to F r {x). 

For r = or r = 1, /x r = and the two independent 
solutions of ea. fT3")> are Fq(x) = 1 + x n and Fi(x) = 1 — 
x N , which correspond to the two degenerate eigenvectors 
|eto) and |<zi). For r = 2 and r = 3 the solution can 
also be found explicitly; the general formula can then be 
extrapolated from these simple cases. We find 



In order to write down the right side of eq.(JHJ we define 
the auxiliary functions 



N-l 



lr { x ) = 



A 



N-l 



im x m 



Til 



h r (x) = Y 



A 



rm 



N — m 



(9) 



It is easy to check that dg r /dx — f r /x and dh r /dx = 
Nh r /x — fr/x. After multiplying Eq.© by x m and sum- 
ming over m, we differentiate both sides with respect to 
x to obtain 



d 
dx 



{1-xf 



fr(x) - A rN - X X 



N 



= - — h r (x). (10) 
X 



with 



F r (x) = (I - x) 1 



r-l 
^ dj"p x 

v=i 



(14) 



_ B{r-l,p) B(N + r-l,p) 

rp 1 ' B{N-l,p) 

Finally, the coefficients of the r-th eigenvector are 
given by 



I d m F r {x) 
to! dx m 



The solution of the differential equation for h r can 
be obtained in terms of its Green function, satisfying 
dG/dx — NG/x — S(x — y). In this case G is given by 
{x/y) N if x > y and zero otherwise. Therefore, 



h r (x) — x N ( a — 



Substituting Eq. lllj l into (|10fl , re-arranging the terms 
and differentiating once again with respect to x, we ob- 
tain 



Since F r is the product of two simple functions, its deriva- 
tive can be calculated explicitly at x = 0. Writing 
F r (x) = N r {x)Q r {x) with 



Q r (x) = (l-x) 1 -' 
we find 



and N r (x) = 1 + 



^ dj"p x 



1 ^ , , d m - p N r 
a rm = — : > B{m,p) 
to! ' 

p=0 



dx m ~P 



dPQ r 



(x=0) 



dxP 



(x=0) 



d_ 

dx 



d ({1-xf 



Sr{x) 



_x N - 1 dx 
Finally, defining 



F r {x) = Q—E£f r ( x ) 
x 

we obtain the differential equation 



pi i _ 1 pi fJ. r F r 



fr{x) 



x (1 - xf 



(12) 



(13) 



Working out the derivatives we find the explicit formula 
valid for r > 2. 



r-l 

E 

p=0 



B{m-p + r-2,r-2)d rp (15) 



for m = 1, 2, . . . , N — 1, with a r g = 1 and a r N = {— l) r . 

From the recursion relations we find that the coeffi- 
cients of the left eigenvectors are given by 



a rm [m{N - m)/N] 



(16) 
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FIG. 1: Transition probabilities for a network with N = 101 
and p = 0.1. The lines correspond to our theoretical cal- 
culation (thick) and to simulations for random (dotted) and 
scale-free (thin) networks, both with an average of 6 connec- 
tions per node. The numerical probabilities were computed 
running the simulations 2 x 10 5 times. The dotted line is 
nearly indistinguishable from the thick line. 



for m — 1, 2, . . . , N — 1, with b r o = b r ^ = 0. Finally, the 
normalization factors T r can also be obtained explicitly: 



rv = 



-! B(N - 



l,r) 



(2r-l) B(N,r) 



(17) 



for r = 2, . . . , N and T = 2, r x = 2N. 

We are finally in a position to state some important 
results concerning the dynamics. The transition proba- 
bility of starting the network at time zero with n = M 
and finding it at a later time t with n = L can be com- 
puted using Eqs.0, 10, lfT5|). iflU and JT7J: 



P(L,t;M,0) = 



N—M 



N 



3M(N-M) , f 
(N+1)N A 2 



JL0 



M 
N 



3M(N-M) ,f 
(JV+l)jV A 2 



Sln + 



6M(N-M) r-i x \n X \\t , v^JV 1 i \i 
N(N 2 -1) \ L ~°L0)[ 1 - — °LN)A2 + <Lr=3 fr rM rL r' 



Since all eigenvalues (except for Ao and Ai) are smaller 
than one, in the limit of long times the transition prob- 
ability is dominated by Ao = Ai = 1 and by the 
largest non-trivial eigenvalue, A2, whose contributions 
we have written down explicitly. For large networks 
we can approximate A| « exp — [2(1 — p)t/N 2 ), so that 
the characteristic duration of the transition process is 
t = N 2 /(l —p), which increases with the square of the 
network size. 



The only two possible asymptotic states are n = 
and n — N, whose transition probabilities from an ini- 
tial state n — M are (N - M) /N and M /N respec- 
tively. A typical transition probability P(L,t;M, 0) for 
L, M ^ 0, N starts at zero if L ^ M, reaches a maximum 
and decreases back to zero. This represents the proba- 
bility that a perturbation initially affecting M nodes will 
lead to a response by L nodes at a time t later. Figure 
1 shows an example for a network with N = 101 nodes 
and p = 0.1. The estimate t 11000 works well for all 
the transition probabilities shown. These results show 
that the theoretical model for fully connected networks 
is an excellent approximation for the average behavior 
of sparsely connected random networks. The theory also 
reproduces qualitatively the behavior of scale-free net- 
works. The deviations from the theory in this case re- 
flect the significant topological differences between the 
two networks. As a final remark we note that the case of 
more internal states per node can in principle be treated 
in a similar fashion. However, although Eq.Q can be 
easily generalized, Eq.@ would be more complicated, 
since it would take more than a single integer n(x) to 
describe the network state. 
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